Bosonic gas as a Galactic Dark Matter Halo 
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We study in detail the properties of gravitationally-bounded multi-state configurations, made of 
spin-zero bosons, in the Newtonian regime. We show that the properties of such configurations, 
in particular their stability, depend upon how the particles are distributed in the different states 
they are composed of. Numerical techniques are used to distinguish between stable and unstable 
solutions, and to determine the final configurations they evolve towards to. Multi-state equilibrium 
configurations can be used as models of galactic halos made of scalar field dark matter, whose 
rotation curves appear more realistic than in the case of single-state configurations. 
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I. INTRODUCTION 

It has been known for a long time that, within the 
context of Einstein's General Relativity, the luminous 
matter content of galaxies cannot explain the so-called 
Rotational Curves (RC) [H-i] , which are still considered 
one of the cornerstone evidence for the existence of non- 
baryonic dark matter. There are many candidates for 
dark matter particles, the most popular ones are known 
as weak interactive massive particles (WIMPS) [4-6]. The 
accepted paradigm that describes the way in which those 
particles form structures is the so called Lambda Cold 
Dark Matter (ACDM) model0-(3- 

An interesting alternative some of us have been work- 
ing on is to consider a (real) scalar field as a dark matter 
candidate, a hypothesis that has been widely explored 
in the specialized literature by many other authors |10l- 
l25j |: see also [IB] for a comprehensive review. In most 
scalar field models, the dark matter particle is an ultra- 
light massive boson, with a Compton wavelength of as- 
trophysical proportions and a very large mean number 
density, so that their collective behavior is well described 
by a classical scalar field. The Scalar Field Dark Matter 
(SFDM) model, as we shall call it in general, offers the 
same results as the ACDM model at large scales, up to 
linear order perturbations 0, [U, [H [H [27|, H| . 

The RC problem has also been addressed using scalar 
fields, see for instance [T^-O EE HI). These works 
considered the dark matter halo as a Newtonian Bose- 
Einstein Condensate (BEC), in which the scalar field dy- 
namics is driven by the so-called Schroedinger-Poisson 
system of equations. However, none of the studies car- 
ried on so far have not shown, undeniably, that these 
scalar field models can account for all features of realis- 
tic galactic halo. 

The modeling of scalar field halos was based on the 
(nodeless) ground state solutions of the SP system, which 



is the only stable solution, and the predicted rotation 
curves are marginally in agreement with the observed 
onespH4l3l I30l - l32j ]. In all cases above, the only stable 
scalar field configuration is that in which boson parti- 
cles are all in the ground state. The ground state is 
the only stable solution of the SP system against grav- 
itational perturbations; other excited configurations are 
intrinsically unstablefT^, [HJ (see alsofH, Hj] and refer- 
ences therein). 

The main purpose in this work is to further explore the 
proposal that was first put forward by Matos & Urena- 
Lopez in Ref. 35 1: that realistic scalar field galaxy halos 
must be comprised of multi-state configurations. As we 
shall show, equilibrium configurations of the SP system 
can be constructed in which many-particle states coexist 
simultaneously, so that the whole system is stable under 
small (radial) perturbations^. Probably not surprisingly, 
we have found that RC could be better fitted by these 
many-particle systems. 

The plan of the paper is as follows. In Sec. [Til we 
present the mathematical theory behind multi-particle 
states. In Sec. IIIII we show that their general properties 
depend upon the distribution of the particles in the differ- 
ent excited states. In Sec. lIVl we give numerical evidence 
that there are stable configurations under small radial 
perturbations, and investigate the late time behavior of 
unstable configurations. The RC curves predicted by sta- 
ble multi-state configurations are calculated in Sec. I VII 
Finally, some conclusions are given in Sec. IVII1 



II. MATHEMATICAL BACKGROUND 

Here we give a brief description of a gravitational 
bounded system of self-gravitating scalar field particles 
following the argumentation in the seminal paper [38j, see 
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1 The relativistic version of the multi-state hypothesis was studied 
recently in|36i, [37ll . in which stability was also confirmed. 
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also[33j, 134 |39H41| . We pay special attention to the key 
aspects needed to build systems that have particles in 
the ground state but also in the excited states (35^. 

We start by assuming a spherically symmetric metric 
of the form 



ds 2 



-a 2 (t, r)dt 2 + a 2 (t, r)dr 2 + r 2 dil . 



(1) 



in units of h = c = 1. The many boson-system is then 
described by a (secondly quantized) real scalar field op- 
erator of the formWa 



X] [hlm$nlm{t, x) + b\ llm <5>* nlm (t , x) 

nlm 



(2) 



where S„; m and b , are usual annihilation and creation 
quantum operators, which obey the commutation rela- 
tions 



bnlm, b^'l'm' 
bnlm 3 bnlm 



= 0. 



(3a) 
(3b) 



The field coefficients $ T iim satisfy the Klein Gordon (KG) 
equation in a curved spacetime 



(□ -,u 2 )$„ /m (i,x) =0, 



(4) 



with □ = (l/y/—g)dfj, (yj—g <9 M ) the covariant 
d'Alembertian operator, and /i is the mass of the 
bosons. The most general solution of Eq. (j4]) is of the 
form 



§ nlm (t , x) = Rni (t, r)Y lm (9,cp), 



(5) 



where we R n i(t, r) is the radial function to be determined 
from the KG equation. The scalar product of the func- 
tions above is defined as 

(^nlm,^n'l'm') = ~i $nlmdnK>i>m' ri? a^y dE , (6) 

JY, 

where 7 is the determinant of the 3-dim metric on the 
spacelike hypersurface E, is a timelike unit vector or- 
thogonal to E, and c£E is the volume element. In our 
case, see Eqs. (JTJ and J5]), Eq. © reads 



($nlm,$n'l'm') = -ifoi'Smm' / R n id t R n ,i,a {t,r)r dr , 

Jv 

(7) 

where we have made use of the orthogonality condition 
of the spherical harmonics Yi m . 

Assuming that there exists a vacuum state defined by 

S nJm |0,0,...,0}=0, V(n,l,m) (8) 
we can construct the orthonormal many-particle states 
\Q) = \N 10Q ,N 2 oo,N 21 - U N 2W ,...), (9) 



JViooWaoo! • ■ ■ 



-|0,0,...0) 



composed of many scalar particles distributed in sets 
of l N n i m particles of mass [i with angular momentum I 
and azimuthal momentum m; the n sub-index labels the 
eigenstates according to their radial function R n i. No- 
tice that many-particle states are constructed from the 
vacuum through the repeated application of the creation 
operators b'. 

On the other hand, the gravitational field is a clas- 
sical field whose dynamics is described by the Einstein 
equations 



G aP = 8irG(Q\ : t a p : \Q) 



(10) 



where the source on the r.h.s. is the expectation value of 
the energy-momentum tensor operator 



T a0 = djzdp® 



d a <$>d a <S> + n 2 ® 2 



(11) 



Notice that we are implicitly assuming the so called nor- 
mal ordering operation, : T a p :, so that 

: b n i m b\ dm := P nl J>nlm , (12a) 
(0,...,0,0|:f a/J : |0,0,...,0>=0, (12b) 

and then the (divergent) vacuum energy density identi- 
cally vanishes. 

The orthogonality of the quantum states ensures that 
the expectation value is given as a superposition of the 
expectation values of the energy-momentum tensor com- 
ponents for each individual state, that is, 

OO 71—1 I 

(13) 



n—1 l—l in— — I 



where we are defining the single states 

I AW) = |0,0,...,A^„ /m ,0,...,0). (14) 
Hence, the Einstein equations (fTU|) read 

(15) 



where 



T, 



af3(nlm) ~ 9 a <& n l m d p<&* nlm — -<? a /9 (9"^nim9 ff $* i m 



Ira) i 



(16) 



and we also normalized the eigenfunctions so that $ — > 
$/V2AW. 

Therefore, in the case when particles populate various 
excited levels, the source of the Einstein equations (fT5|) is 
equivalent to the energy momentum tensor of many (in- 
dependent) classical complex scalar fields $ n im(t, x) min- 
imally coupled to gravity. Each one of such scalar fields 
accounts for only one of the excited single states (IT41 . and 
its dynamics is given by its own KG equation (Q} [3a, |42J . 
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Finally, we consider the Newtonian limit of the coupled 
Einstein-Klein-Gordon (EKG) equations (gj and ([15]). 
which results in the so-called Schrodinger-Poisson (SP) 
system (43j 



v 2 [/ = V|f„ !m | 2 



nlra 

1 

~2 

where \& n 7 m is related to <£>, 



Ira 



V 2 ^ 1 

v 1 nlm 



by 



(17a) 
(17b) 

(18) 

Then, the Newtonian version of the EKG equations de- 
scribes the dynamics of non-relativistic wave functions 
which are coupled among themselves through the New- 
tonian gravitational potential U . 

Once in the non-relativistic regime, we can define phys- 
ical quantities like the kinetic K and gravitational W en- 
ergies, and the total number of particles Af. These quan- 
tities can be explicitly given in terms of the Newtonian 
fields as 1331 



K 



n.l.m 



V nlm v n 



Ira H~ 



W = J U\^nlm 



2 dv. 



^= E 



l dv. 



)4i9a) 

(19b) 
(19c) 



These expressions will be useful later to monitor the nu- 
merical evolution of multi-state configurations studied in 
Sec. El 



III. MIXED NEWTONIAN STATES 





<Mr) 
4>,(r) 



«|> 2 (r) 
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FIG. 1. Radial functions 0i and 02 for the mixed configura- 
tion |JVi, 1.1), and </>i, 02 and 03 for the mixed configuration 
|iVi, 0.96, 0.91) see also the first and third entries in Table U 
and text below for details about its construction. 



In this section we construct solutions of the SP sys- 
tem (TTT1) when JV bosons are allowed to occupy 1 differ- 
ent levels, all of which, for simplicity in the discussion, 
will have zero angular momentum (I = 0, m = 0). Hence, 
the states are of the form \Q) = \N X , N 2 , N3, N x ). As- 
suming spherical symmetry, we then have 

1 8 2 (r 2 U) 1 



dr 2 



E^ 



1 a 2 (r 2 *, 



U9 n , 



(20a) 



l,42Db) 



dt 2r 2 dr 2 

First of all, we look for stationary equilibrium config- 
urations in the form 



* n = e-^Vr»(r-), 
for which the system (|20|) becomes 

1 d 2 (r 2 U) 



r 2 dr 2 

1 rf 2 (r 2 0») 
r 2 dr 2 



= Ei^i 2 ' 

n=l 

2{U - 7n )0„ , 



(21) 

(22a) 
1,..,Z (22b) 



Eqs. (|22p will be solved under the following boundary 
conditions. In order to obtain regular solutions at the ori- 
gin, we demand that the spatial derivatives of all <p n and 
U must be zero at this point, but we arbitrarily prescribe 
the central values of all 0„, as these values are the free 
parameters of the solutions. Also, we impose n (?*) 
and U(r) = —Af/r as r — > 00 because we are looking for 
bounded configurations. 

With these boundary conditions the system becomes 
an eigenvalue problem. Given the central values {0„(O)} 
there are unique values {jn} and U(0) for which the 
boundary conditions are satisfied. The numerical solu- 
tions are then found by using a shooting method to inte- 
grate Eqs. (|22p from r = up to the numerical boundary 
f = f m ax, with {"f n } and U(0) playing the role of shoot- 
ing parameters. 

Because the equations are integrated in a finite nu- 
merical domain, it is more convenient to introduce more 
detailed specifications to the boundary conditions at 
r — > 00. We use the asymptotic behavior of the solu- 
tions: the gravitation potential goes like U(r — > 00) = 
—Af/r, and the radial functions behave as 4> n {r — > 00) ~ 
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State 


040) 


7n 


K 


W 


N 


|^i,l.l) 
|iVi,1.6) 
\Ni, 0.96, 0.91) 


2 (O) = 0.756 
02 (0) = 0.934 
2 (O) = 0.710, 3 (O) = 0.543 


7i = -1.033, 72 = -0.574 
7i = -1.163, 72 = -0.677 
71 = -1.185, 72 = -0.712, 73 = -0.471 


1.846 
2.272 
2.407 


-3.691 
-4.544 
-4.819 


3.493 
3.945 
4.520 



TABLE I. Central values of excited states 0n(O), eigenvalues y n , kinetic K and gravitational W energies, and the total number 
of particles J\[ of three mixed states, all of them with 0i(O) = 1.0. The labeling of the mixed states is in terms of the r\ 
parameters defined in Eq. (|25[) . See text bellow for details. 



exp(— yj— 2j n r). Then, more suitable boundary condi- 
tions are 

U{r max ) + T max U max = , (23a) 

(t)'n( r max) + V^l?i 4>n (r m ax ) = . (23b) 

The shooting procedure is then used for different values 
of r max . As r max is increased the shooting parameters 
converge, and we choose as solutions those which satisfy 
the boundary conditions (|23p within a prescribed toler- 
ance. 

Before we show the numerical results we make some re- 
marks about notation. In the case when all the particles 
are in the one same level, the system (l22l consists of only 
one Schrodinger equation and only one wave function in 
the source term of the Poisson equation. Such system 
has been widely studied in j32| . where several solutions 
{4>ni 7n> U } have been calculated. 

The main difference between the functions n is that 
they have (n — l)-nodes in their radial profile. The state 
corresponding to the zero-node function 0i has the lowest 
(negative) total energy E — K + W, the lowest energy 
eigenvalue 71, and is correspondingly the less massive; 
this state is called the ground state. The other solutions 
with nodes are more massive and have larger energy val- 
ues than the ground state; they are called excited states. 

In this spirit, we shall call mixed states to those config- 
urations in which single states <p n are present simultane- 
ously. For instance, the simplest mixed state configura- 
tion consists of the ground state 0i and the first excited 
state 02, and then it is characterized by quantum state 
\N\, JV2). In the next section we construct some of these 
configurations. 

A. Ground-first mixed states \Ni,N2) 

Let us start with the simplest mixed state, that com- 
posed of the ground state 0i plus the first excited state 
02- At this point, it is not necessary to solve Eqs. (|22|) 
for each possible pair (0i(O), 02(0)) in order to find the 
complete space of solutions. 

Instead, it is worth to use the scaling symmetry that 
the complete SP system obeys (32j; for the particular case 
here, it reads 

{0i,02,7i>72,*/,r} -> (24) 



where A is an arbitrary parameter. 

This means that once we have found a solution to 
the SP system for given values of (0i(O), 02(0)), there 
is a complete set solutions each of which are related to 
each other just by the scaling transformation (IM|) . This 
set we will call it a family of solutions. Different fami- 
lies are then found by taking different central values of 
(0i(O),«fc(O)). 

The central values 0i(O),02(O) of the solutions in a 
family, as a consequence of the scaling relation (PMl) , will 
be located along the straight line defined by the origin 
and the given point (0i(O), 02(0)) on the plane 0i,02- 
Therefore, once all the solutions with the same value 
0i (0) and different 02(0) are known, and their respec- 
tive families have been calculated, the complete space of 
solutions can be constructed as the collection of all fam- 
ilies of solutions. For simplicity in the notation, we will 
drop the caret symbol from the field quantities, in the un- 
derstanding that we are dealing with scaled quantities. 

Useful quantities to characterize mixed states are the 
ratios of the number of particles in different states with 
respect to the ground state; we define these ratios as 

Vn = A/„/.A/i , (25) 

where by definition 771 = 1, and the total number of par- 
ticles is Af = A/"i(l + r/ 2 + r/3 + • • • ). Notice that r] n are 
invariant quantities under the scaling relationship (|24[) . 

We calculated solutions with 0i(O) = 1 for different 
values of 02(0), and we observed that 772 increases mono- 
tonically as 02(0) grows, starting at the value 772 = 
when 02(0) = (no particles in the excited state, N 2 = 
0). This behavior of 772 allows us to choose it as a free 
parameter instead of 02(0), and then (0i (0)5/72) will be 
the representative parameters in the construction of the 
solutions. This is a convenient option because, as we will 
see in the next section, there is evidence that the sta- 
bility under radial perturbations of a ground-first state 
depends mainly on the values of 772 . 

In order to fix the value of 772 , it is necessary to add the 
differential expressions for the number of particles Afj 

^ = <t>V, n = !,..,!, (26) 

in the system (j2"2")l . The boundary conditions for these 
equations are given at the origin N n (0) = 0, and the 
desired value of 02 is imposed at r max . We then solve 
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the system of equations (1221) and (l26l) using a shooting 
method; this time, however, the shooting parameters are 
%(0), 71, 72 and t/(0). 

Once again, the complete space of configurations is 
formed by the collection of families of solutions, each one 
characterized by the same value of <j)\ (0) and different 772 . 
As a consequence of (fMf , the physical quantities of the 
solutions for each family will be related by 

{N X ,N 2 ,K,W}^ (27) 
{\N U \N 2 ,\ 3 K,\ 3 W} . 

In the top panel of fig. Q] we show typical radial func- 
tions of a (mixed) ground-first state, the zero-node ra- 
dial function corresponds to the ground state whereas the 
one-node radial function corresponds to the first excited 
state. The mixed state was constructed with 0i(O) = 1.0 
and 772 = 1.1, and so we labeled it as \N\, 1.1). 

In the first row of Table HI we show the scalar field 
central value for the excited state 02 (0), the frequencies 
for both states, 71 and 72, as well as the kinetic and the 
gravitational energies, and the total number of particles 
of the system. The same quantities are shown in the 
second row of the same table for the mixed state |A/i, 1.6), 
which was constructed with </>i(0) = 1.0 and 7/2 = 1.6. 

Following a similar procedure to the explained above 
we also constructed systems where particles are coexist- 
ing in the ground state and in the first and second excited 
states. In the bottom panel of Fig. [TJ we show the radial 
functions for one of those systems with 0i(O) = 1.0, and 
m = 7/3 = 1.0; we called it \Af u 1.0, 1.0). 

In Table U we show its energy eigenvalues 71,2,3 and 
other physical quantities. We were able to verify that, 
within the limitations imposed by the numerical error, 
all configurations we could find satisfy the virialization 
condition 2K + W = 0. 

We constructed several ground-first states for different 
values of ^i(0) and r/ 2 . In Fig.[3l we plot the total energy 
of each system in terms of the number of particles in 
the ground state and first excited states, Ni,N 2 . It is 
possible to notice that the total energy for each system 
is negative, which implies that they are gravitationally 
bounded objects. 

In Fig. HJ the energy eigenfrequencies 71 and 72 are 
shown; for all configurations 71 < 72- Finally, in Fig. [SJ 
the kinetic energy of each state, K\ and K 2 , is shown. In 
contrast with the behavior of the eigenfrequencies, there 
are configurations for which K\ < K 2 if 772 > 1, whereas 
Ki > K 2 if 772 < 1. 



B. Mixed states as a generalization of single states 

We want to show now that mixed states are gener- 
alizations of single states. As mentioned before, if the 
number of particles M in an single state is fixed, what- 
ever the ground or any of the excited states, the total 
energy E takes a fixed value because of the eigenvalue 




/ • - - E (ground state) 

• ; — E (first state) 

-0.6 '- : - - - E (ground-first state) 

I : ■ - ■ y 2 (first state) 

.... y, (ground-first state) 

-0.9 -i 



-1.2 -j 

J 1 I 1 I 1 L 

10 20 30 



FIG. 2. Comparison of the total energy E = K + W of a 
ground-first configuration with a fixed number of particles 
M = Mi + M2 = 2.0622, with two single-state configurations 
with the same number of total particles: one ground single- 
state, and one first excited single-state. The total energy 
depends upon the fraction r/2 = M2/M1, see Eq. (f25)l . If 
Mi 3> M2, E goes to the energy value of the ground single- 
state; if Mi <C A2, then E goes to the energy value of the 
first excited single-state. We make a similar comparison with 
72. If Mi 2S> M2, 72 takes very large negative values, whereas 
if Mi <C M2 , then 72 approaches the value corresponding to a 
first excited single-state. 



problem; this value however increases for larger number 
of nodes in the radial profile. 

This behavior of E changes radically if the same num- 
ber of particles M are distributed in a mixed configura- 
tion. Now, the total energy of the system E can take 
values from a continuum interval depending on how the 
particles populate the states of the mixed state. 

In particular, for a ground-first configuration the total 
energy is a monotonically increasing function of 772- If 
77 ~ 0, E — > £;( slnglo ) ; where E\ is the total energy of 
a single ground state composed of M particles. On the 
other hand, if the particles are moved from the ground 
to the first excited state, the total energy of the mixed 
state takes continuum values, and E — > E^ 111 ^, where 
^(single) j g .j-Q-j-gj ener gy f a single first state, as the 
ground state is depopulated (in other words, 772 — > 00). 

In Fig. [5] we show this behavior of the total energy 
for a system of with a fixed total number of particles 
M = 2.0622. Notice that the values are within the range 

^(single) < E < jingle) ^ ^ q{ 

was found for all other characteristic quantities, like the 
eigenfrequencies 7„, see also Fig. [2j 

From this perspective, we can say that, given a system 
of M particles in a single-state configuration, it is possible 
to change the properties of the quantities attached to 
it by moving particles away to populate other excited 
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Total Energy 




3 0.0 



FIG. 3. Total energy of ground-first configurations in terms 
of A/i and A2. All of them are negative, which implies that 
the systems are gravitationally bounded. 



Energy Eigenvalues 




3 0.0 



FIG. 4. Frequency eigenvalues 71 and 72 of ground-first con- 
figurations in terms of A/"i and A/i- For all the configurations 
it is satisfied that I72 < |7i|. 

states. 



IV. STABILITY OF MIXED STATES UNDER 
SMALL PERTURBATIONS 

In order to promote the existence of mixed states be- 
yond the mathematical context, it is first necessary to 
prove their stability. It is known that when all the par- 
ticles are in the ground or in an excited state, the con- 
figuration is stable under small radial perturbations that 
strictly conserve the number of particles, SN = (l7ll3^. 

Instead, if the system is considered open, so that SN — 
is not demanded, an excited state emits particles, loses 
its nodes, and eventually settles down on to a ground 
statejl^Hl]- On the contrary, grou nd states can tolerate 
perturbations for which SN ^ 0|44j, and it is in this sense 
that ground states are said to be stable, whereas excited 
states are unstable. 

It is then expected that the excited state of a ground- 
first configuration should remain unstable if there are 
very few particles in the ground state. However, the re- 
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FIG. 5. Kinetic energy for each separate state K\, K2, in a 
ground-first configuration in terms of Mi and A/2- The sys- 
tems with A"i > K2 correspond to 772 < 1 (A/i > A/2), whereas 
those with K\ < K2 correspond to 772 > 1 (A/i < A/i). 

suits in the previous section show that the properties of 
excited states change depending on how many particles 
populate the ground state. 

This is actually the case of stability. We found that, 
by adding particles to the ground state, it is possible 
to construct ground-first configurations for which, under 
open conditions, the mixed state is stable under radial 
perturbations. As particles do not interact directly one 
with each other, the change in the stability of the excited 
state is produced just by the gravitational interaction 
with the ground state. 

Stability studies of single configurations have been 
done with perturbation theory and full numerical evo- 
lutions. In principle, the stability of mixed states can 
be also studied perturbatively (a work that is worth a 
separate manuscript). Instead, we have chosen to evolve 
numerically mixed states not only to determine their sta- 
bility, but also their late time behavior. We have done 
this for several ground-first configurations with different 
values of 772 • 



A. Numerical perturbation 

The evolution of the mixed states is done by solving the 
discretized version of the time dependent SP system (|20p , 
taking as the initial data the functions of the different un- 
perturbed states ^n, see Eq. (f21"T) . constructed in the pre- 
vious section. We consider no perturbations other than 
those introduced by the finite differencing error in the 
numerical integration. 

The procedure followed in the construction of the so- 
lutions to Eqs. (|20|) can be summarized as follows. 

1. The numerical grid is populated with the initial 
data 

tt*(*o =0,fcAr) = </> n (kAr), 
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where to is an arbitrary initial time that we choose 
to be zero, Ar is the spatial resolution of the grid, 
the super-index N labels the discretized numerical 
solution, k is an integer that labels the grid and n 
runs from 1 . . .X, X the number of populated states. 

2. The gravitational potential U N (0, kAr) is obtained 
by introducing ^ in Eq. ([20]). 

3. Using the obtained gravitational potential U N 
in Eqs. (|20|) . each populated state of the sys- 
tem is leaped forward in time a step At getting 
^ (At, kAr). 

4. Repeating j — 1 times the steps 2 and 3 above we 
obtain V%(jAt, kAr) and U N ((j - l)At, kAr). 

Because the SP system is discretized in order to solve 
it numerically, it is expected that the numerical solution 
differs from the equilibrium value ^ n in Eq. (|2"TT) by 
a numerical error AvE'n, i.e., 

*%(jAt, kAr) = e - l7 " 3A Vn(fcAr) + AV n (jAt, kAr) . 

We then say that the system is perturbed without con- 
sidering an explicit perturbation, because the numerical 
error A^J n , that comes from the discretization, is consid- 
ered as the perturbation itself; in fact, behaves like 
a perturbed ^ n during the evolution. 

We have to stress that during the numerical evolution 
we allow the system to eject particles. With the imple- 
mented open boundary conditions in the code, the condi- 
tion for which the number of particles has to be preserved 
is not maintained. 



B. Perturbing the ground state 

In Ref. [HI the discretization error was used to perturb 
a single ground state in order to study its stability and it 
was shown that the numerical evolution reproduces the 
results obtained with perturbation theory. Here we give 
a brief account of this result. 

Using perturbation theory to first order, when a single 
ground state function = 0i(r)e -47lt is perturbed with 
a small radial perturbation S^i, keeping the number of 
particles constant, a new oscillation mode a appears in 
the perturbed system ^ pert = ^\ + 8^\. It is found 
that <5M/i is regular, spatially localized, and has an har- 
monic time dependence (such time-dependence involves 
not only the new oscillation mode a, but also the eigen- 
value 71). The system is then said to be stable under 
small radial perturbations. 

A quantity that gives relevant information about the 
new oscillation mode is the perturbed density. The 
ground state non-perturbed density p\ = \^\\ 2 is time 
independent, whereas the perturbed one p per t = \^pert\ 2 
has an harmonic time dependence and oscillates with an 



angular frequency 2na. To first order in S^i, the per- 
turbed density is 

Ppert = |*i + S^\ 2 = Pi(r) + Sp(r) cos(27rcrf) . (29) 

These perturbed ground state properties are also found 
when the ground state is evolved numerically without 
considering any explicit perturbation, except for that in- 
herent to the discretization. Using as initial data the time 
independent function cf>i of a ground state, its temporal 
behavior is obtained as explained in the previous subsec- 
tion. A quantity that is monitored throughout the evolu- 
tion is the function Re[4'^ v (0, t)}, which behaves harmon- 
ically in time as expected. Its Fourier Transform (FT) 
shows a main harmonic mode that matches the eigenfre- 
quency 71 of the unperturbed ground state. 

The density p?(t,Q) = |*f (i, 0) | 2 also shows an har- 
monic behavior in agreement with the perturbative re- 
sult (|29p . It oscillates around the value of the non- 
perturbed density pi, and from its FT we observe that 
its main oscillation mode coincides with a. Finally, even 
with the open boundary conditions in the numerical code, 
the number of particles does not change, which is consis- 
tent with the results of perturbation theory. 

C. Perturbing stable mixed states 

In order to study the stability of a mixed state, we 
evolve it following the steps described in IIV Al The 
general idea is to compare the behavior of Re[*^(t,0)], 
Pn (t, 0), and the number of particles Af„ (t) for each oc- 
cupied state with the behavior of the corresponding quan- 
tities for a single ground state. We can conclude that 
there is evidence of the stability of the mixed configura- 
tion if those behaviors are similar. 

In fig. [6l we show the numerical evolution of 
Re[^^(t, 0)] for a ground-first configuration with 772 — 
1.1, i.e., a |A/i,l.l) state. The wave functions behave 
harmonically, and their FTs, presented in Fig. [7J show 
that the main harmonic mode of each corresponds 
to the angular frequency j n of the unperturbed *5f n . 

Furthermore, in Fig. [5] we show the numerical values of 
the central densities p% = \^f^ 0)| 2 , an d it is clear that 
they oscillate closely around the value of the unperturbed 
densities (which are formally time- independent), p\ — 
|</>i| 2 = 1.0 and p 2 = \fo\ 2 = 0.572. 

The bounded oscillations suggest that the numerical 
perturbations of each ty n are spatially localized and have 
an harmonic time dependence. In the bottom panel of 
Fig. [51 we show the central value of the total density 
p N = pi + P2 , with its corresponding FT in Fig. [9] As 
in the case of a single ground state, for which its quasinor- 
mal mode is the mean harmonic mode of the perturbed 
density, we expect that the mean harmonic mode shown 
in this figure corresponds to the characteristic oscillation 
mode of the perturbed state \Afi, 1.1). 

In Fig. 1101 we verify that the number of particles for 
each occupied state is conserved separately, so that the 
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FIG. 6. Evolution of Re[*f (t,0)] (top) and Re[*f (t, 0)] (bot- 
tom) for the mixed state | iVi , 772 ) - An harmonic behavior is 
observed, and the oscillation modes can be read from the cor- 
responding FFT shown in Fig. [7] 



same happens for the total number of particles. In the 
bottom panel of Fig. [101 the relation 2K + W for the 
system is presented, which in turn shows that the config- 
uration remains virialized during the evolution. 



D. Perturbing unstable mixed states 

We performed the evolution for several ground-first 
states with <px(0) = 1.0 and different values of 772- A 
general result is that if r\2 < 1.1 the mixed state behaves 
similarly as the state \N\, 1.1). That is, each state of the 
configuration evolves harmonically with the main angular 
frequency j n of the unperturbed wave function the 
oscillations of the central densities p 1 ^ are bounded, and 
the number of particles in each state is conserved. We 
take for granted that these characteristics are evidence 



FIG. 7. (Top) FFT of Re[*f (t,0)] for the |iVi, 772) mixed 
state. The main mode at f — 0.164 corresponds to the angu- 
lar frequency 2nf = 1.033 which coincides, in good approxi- 
mation, with the value of 71 of the unperturbed wave function 
*i, see TableU (Bottom) FFT of Re[*f (t, 0)] for the \Ni,ip) 
mixed state; the main mode at f — 0.091 corresponds to the 
angular frequency 2nf = 0.574, which also coincides with the 
value of 72 of the unperturbed wave function $2, see Tabled] 

of the stability of these ground-first configurations. 

On the other hand, the evolution of configurations with 
r/2 > 1 . 1 differ considerably. In the top panel of Fig. fTS] we 
show the early evolution of the central densities p% for a 
ground-first configuration with </>i(0) = 1.0 and 772 = 1-6, 
which is the state \N-y, 1.6). 

At the beginning, they just oscillate around the val- 
ues of the unperturbed initial state, pi(0) = 1.0 and 
P2(0) = 0.873. However, after a while the amplitudes 
of the oscillations start to grow. In fact there is a pe- 
riod of time in which it is possible to fit each p 1 ^ with a 
function of the form 

p n ~ aie Q2 * cos(a 3 i + 04) , (30) 
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FIG. 8. (Top) The perturbed central densities p™ = 
0)| 2 for the [iSTj , 772) configuration are shown. They os- 
cillate around the constant density values of the unperturbed 
states, corresponding to p\ — 1.0 and p2 = 0.485. Such os- 
cillatory behavior is the expected one for stable states (see 
Sec. IIV Bl for details). (Bottom) Overlap between p N (Ar) 
and 4(p JV (Ar/2) - p N (Ar/4)) + p N (Ar/2). This plot shows 
that the numerical evolution is second order convergent. 



where the a's are constants. This behavior suggest that 
the perturbations of the states \& n , besides their har- 
monic time dependence, are exponentially unstable since 
a 2 > in general. 

From Fig. [TJJ it is possible to infer that the exponen- 
tial growth of the density amplitude of the excited state, 
which starts almost from the very beginning of the evo- 
lution, acts as a trigger for the instability of the com- 
plete configuration. For completeness, we show in the top 
panel of Fig. [TTJthe convergence of Pi 2 (t, 0) for different 
spatial resolutions. We can conclude that the exponen- 
tial growth of perturbations is a physical characteristic 
and not a spurious numerical result. 



FIG. 9. FFT of the numerically-perturbed total density 
p N (t, 0). The main mode appears at / = 0.065, which should 
correspond to the quasinormal mode of oscillation of the per- 
turbed configuration. 



Moreover, the number of particles in each state of 
the configuration \Ni,1.6) is not conserved. In Fig. [14] 
we show that both the ground and the first excited 
states lose particles. Even more, the behavior of each 
Re[5'^ r (t, 0)] also shows that their angular frequency 
change with time, and that only at the early stages of the 
evolution the frequencies coincide with the unperturbed 
values 71 = —1.163 and 72 = —0.677, see the top panel 
of Fig. [13] (the way in which the plotted quantities are 
computed is explained in the next section). Therefore, 
the configuration |j\Ti, 1.6) is not stable. 

We performed the evolution for several configurations 
with different values of 772 > 1.2, and we noticed that the 
speed of growth of the perturbation amplitudes increases 
for larger values of 772 . In order to quantify the speed 
of growth in terms of 772, during the early times of the 
evolution we fit the central density of the excited state 
P2 with a function of the form (f30|) . 

The values of coefficient a<i for the exponential coef- 
ficient in terms of 772 are plotted in Fig. [TUJ A linear 
extrapolation indicates that a2 — > as 772 — > 1.13. As we 
already found that configurations with 772 < 1.1 do not 
exhibit exponential growth, then we take 7y* rcsh ~ 1.13 
as a threshold value that separates stable and unstable 
configurations. 

This result has been obtained from configurations with 
4>i (0) = 1.0 and different values of 772. Since mixed con- 
figurations obey the scaling symmetry (|24p , for which 772 
is an invariant quantity, then this threshold value should 
hold for all configurations. 
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FIG. 10. (Left) Number of particles for each state of the 
mixed configuration \Ni,f]2) is shown. It is observed that the 
number of particles is conserved all along the numerical evo- 
lution. (Right) The value 2K + W is shown. Since it oscillates 
around zero, it is inferred that the mixed configuration is a 
virialized system 



V. LATE TIME BEHAVIOR OF UNSTABLE 
MIXED STATES 

In this section we present evidence that shows that 
unstable ground-first configurations, when perturbed, 
evolve towards a stable configuration, in a process that 
parallels the evolution of single unstable states into the 
ground one, see^|. 

During the stabilization process, we will see that the 
excited state loses particles and the node in its radial 
profile disappears. Likewise, the ground state loses par- 
ticles as well, but a node appears in its radial profile. In 
any case, it will be possible to identify the final stable 
configuration the system is settle down onto. 

A useful procedure is to follow the evolution of the 



FIG. 11. (Top) The perturbed central density p2 of a unsta- 
ble configuration with rj = 1.6 is shown. The overlap between 
p${Ar) - p% (Ar/4) and 4(pf (Ar/2) - pf (Ar/4)) show that 
p2 converges at second order. (Bottom) Coefficient 0,2 of the 
exponential function (|30[) used to fit the growth in the os- 
cillations of the P2 for different unstable configurations are 
shown. It is null for the threshold value of 772 — 1.13, see text 
for details. 



effective eigenfrequency defined by 



,N _ 



LH*f \ 2 dv 



( 31 ) 

which is a generalization of the eigenfrequency obtained 
from the Schroedinger equation ([20]) in the case of sta- 
tionary solutions of the form (|2Tj) . i.e., 



7i 



<j>iV 2 <j)idv + / U\$dv 



(32) 



The precise identification of the final configuration is 
a tricky task, as the ground and the excited state inter- 
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change the node in their radial profile during the evo- 
lution. We have found that the labeling of the states, 
whether ground or first excited state, depends basically 
in the relative values of 7^ and 7^. As in the case of 
stationary solutions, we shall call ground state the state 
with the lowest value of 7 . 

Now, we are going to describe the main stages in the 
evolution of the ground-first configuration \Ni,l.6) in 
terms of pf (t,0), tf(t), and Aff(t). At the beginning 
of the evolution, pf (t, 0) oscillate with small amplitude 
around the unperturbed values p\ — 1.0 and P2 — 0.873, 
see the top panel in Fig. (fT2|). The same behavior is 
found for 7 l JV (i), which oscillate around 71 = —1.163 and 
72 = —0.677, see the top panel in Fig. [13] 

During this stage 7^ < 7^, we can clearly see that 
the radial profile of Re[$^(t, 0)] shows no nodes, whereas 
Re[<I'^ v (£, 0)] has one node. The number of particles 
in each state remains constant, see Fig. [14] 

Later on, the amplitude of each pf(t,0) starts to 
grow exponentially as discussed before, see for instance 
Fig. [12] Meanwhile, 7^ oscillate with bigger amplitudes 
and move away from their initial values, in such a way 
that decreases and 7^ increases, see Fig. ([LI]) . At 
the same time, the number of nodes in the radial profiles 
Re[\l/f (i, 0)] changes quickly and the number of particles 
of both states starts to decrease, see Fig. ([14")) . 

Then the system enters in a stage in which it experi- 
ences the most violent changes. The amplitudes of the os- 
cillations of pf(t, 0) and 7f (i) becomes larger (Figs, (TT2)) 
and ([T3]) h and the values of 7/ v (t) intersect and continue 
moving away from each other. The number of nodes in 
the radial profiles continues changing, and the number of 
particles of each state decreases quickly, see Fig. ([T4"|) . 

Finally, the system relaxes and the values of pf(t,0) 
and 7i V (t) start to converge and oscillate around a fixed 
value. It can be noticed that 7^ > y% > an d that the 
radial profiles have interchanged nodes. Also the number 
of particles in each state stabilizes around fixed values. 
All this description strongly suggests that the system is 
reaching a stable configuration. 

In order to verify that the final configuration really 
corresponds to a stable one, we read from the numerical 
results a set of parameters that can help us to construct 
an equilibrium configuration as described in Sec. IHI Al 
For the particular case studied here, we find that = 
0.617 (recall that the excited state became the ground 
one), and r/ 2 =Aff/Aff = 0.544. 

We then construct an equilibrium configuration with 
the aforementioned values as input parameters, compute 
all relevant quantities, and compare them with their (nu- 
merically) evolved counterparts. The resulting values 
are shown and compared in Table [TT] Based on the 
fact that those values coincide in good approximation, 
we can affirm that the unstable configuration |Ai,1.6), 
evolves towards a stable configuration with 772 = 0.544 
and 0i (0) 1.557. 

Long evolutions for several configurations for which ini- 
tially 7/2 > 1.1 and 4>\ (0) = 1.0, show that their behavior 



is very much like the one observed in the unstable con- 
figuration with 77 = 1.6. In all cases, they evolve towards 
a stable configuration, as exemplified in Fig. (|15p . 

Here we have defined another scale-invariant quantity 
for the ratio between the eigenfrequencies of the equilib- 
rium configurations, 

T„ = 7n/7i > (33) 

such that again T\ = 1. We can plot the resulting 
values of ground-first stationary solutions on the plane 
(772^2), which are then represented by the solid line in 
Fig. 1151 We notice that equilibrium configurations with 
initial 772 > 1.1 evolve towards an equilibrium configura- 
tion with 772 < 1.1. 
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FIG. 12. Evolution of pf r (t, 0) for a unstable ground-first con- 
figuration with rj = 1.6. (Top) Early time behavior; (Bottom) 
late time behavior. The system eventually settles down onto 
a stable, stationary, configuration. 



VI. ROTATION CURVES IN SCALAR FIELD 
GALAXY HALOS 

We mentioned before that our main motivation for 
studying mixed configurations was the possibility that 
a scalar field would be the dark matter in galaxies. For 
that, we wanted to explore the capabilities of a scalar 
field to form realistic galaxy halos. 

A crude estimation of rotation curves in mixed states 
was first presented in Ref. (35[ , under the (then untested) 
assumption that they were stable, a feature we can now 
consider firmly confirmed by the results of the present 
work. 

For the mixed states studied in the previous sec- 
tions, we can calculate the velocity of test particles mov- 
ing along circular orbits in the gravitational potential 
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Mo) 


02(0) 


7i 


72 


U(0) 


Equilibrium 


1.557 


0.795 


-1.355 


-0.692 


-2.461 


Final 


1.563 


0.786 


-1.355 


-0.690 


-2.461 



TABLE II. Central values of the excited states <f> n (0), eigenvalues 7 n , kinetic K and gravitational W energies, and the total 
number of particles AT of the mixed states. Because of the similarity between the final values of the evolved configuration shown 
in Figs. fT2l IT5l it can be concluded that the final state certainly corresponds to an equilibrium configuration. 




FIG. 13. Evolution of for a unstable ground-first configura 
tion with r\ = 1.6. (Top) Early time behavior; (Bottom) late 
time behavior, see text below. 

sourced by the mixed states configurations via the New- 
tonian formula 

v(r) = y/AT(t,r)/r, (34) 

where Af(t, r) is the total number of particles inside the 
radius r obtained from the numerical equilibrium config- 
urations of the SP system (|2"0)) . 

The results are shown in Fig. [T5] We can see a notice- 
able improvement in the flatness of the rotation curve at 
large radii as long as more excited states arc taken into 
account. Some comments are in turn. 

The circular velocity of mixed configurations shows 
some flat profile at intermediate radii, whereas the typi- 
cal Keplerian tail shows up at large radii, indicating that 
at the end we are dealing with a localized object of fi- 
nite size. We also note the existence of some ripples in 
the velocity profile, which are a consequence of the nodes 
present in the radial functions of excited states fa. 

The height and position of the first peak in the velocity 
profile are approximately set by the ground-state wave 
function fa, whereas the total size of the objects is fixed 
by the profile of the most excited state. 

At this point, we cannot say that mixed states are 
already strong candidates to explain galaxy halos, but for 



FIG. 14. Evolution of Ni and N2 for a unstable ground-first 
configuration with 77 = 1.6. 



that a more complete study would be necessary, like the 
one pursued in |10Hl3l |. where also baryons were included 
in the equations of motion. 

We can though provide some clues about the possible 
physical features of a scalar field galaxy halo. In most 
scalar field dark matter models, the scalar field mass is 
usually very light, around fi ~ 10~ 23 eV; for such a small 
value, it was possible to consider that single state config- 
uration, whether the ground or any of the excited ones, 
accounted for the complete halo configuration. This was 
a central assumption in most of previous references about 
scalar field dark matter models{lE El EE HI S3 • 

The Compton wavelength of the scalar field is very 
large for usual standards, because Ac = — 10 pc. In 
the case of a Newtonian configuration, the size of the 
bounded object scales like R = r/X 2 , where A is the 
scaling parameter in Eq. (|24p . which in turn is related 
to the central field value of the configuration through 
A 2 = <M0). 

Following previous works, the scaling parameter is es- 
timated to be A ~ 10 3 [32], |45j]. This implies that a stable 
single-state equilibrium configuration can model a galaxy 
halo of a size of around 5 — 7 kpc. However, larger galaxy 
halos, as is the typical case, are out of the capabilities of 
single state configurations. This limitation can be no- 
ticed, for instance, in the fits done in[13]. 
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FIG. 15. Fate of unstable configurations with r\ — 1.4,7) = 1.6 
and r\ — 1.8, in terms of parameter F, see Eq. (I33p . Notice 
that the configuration oscillates around the line represent- 
ing stationary equilibrium configurations on the stable branch 
772 < 1.13, see Sec. UVDl 



Mixed states can alleviate this limitation. As the size 
of the configuration is determined by the most excited 
state, the scalar field halos can be as large as necessary, 
and this helps to fit better the RC in real halos. 

Moreover, as already noticed in Ref.j35j], mixed states 
provides us with more free parameters to play with. 
The extra parameters are the occupation numbers of the 
mixed state, namely Ni, N2, Ns, etc. Except for the limi- 
tations imposed by stability, these values would be deter- 
mined initially, before the collapse of the scalar field con- 
figuration, by the local environment a scalar halo could 
be subjected to during the cosmological evolution. 



FIG. 16. The rotation curve v(r), see Eq. (|34[) . for the single 
ground state \Ni — 2.0622), and the mixed configurations 
\Ni, 1.1) and \Ni, 0.96, 0.91), see also TableQ] Notice that the 
flatness of the curve is improved if more excited states are 
taken into account. 



VII. CONCLUSIONS 



We have shown, for the first time, the existence of 
stable many-particle states made of scalar field in the 
Newtonian regime; such states are the generalization of 
the well-known boson stars that have been exhaustively 
studied in the literature. The possibility of mixed states 
was already suggested in the seminal paper of Ruffini 
& Bonazolla about boson stars, but their existence and 
properties had not been studied before 

Detailed instructions were provided for their construc- 
tion and classification, but more importantly is that we 
established simple and sufficient criteria to determine 
their stability properties. Even though Newtonian con- 
figurations obey a scaling relationship, it was possible to 
define invariant parameters that allows us to follow their 
evolution and determine their final fate. 

Some remarks were given regarding the importance 
these results may have for scalar field dark matter mod- 
els in describing the properties of galaxy halos. However, 
more work is needed to have a complete picture of all 
possibilities offered by scalar fields and their gravitation- 
ally bounded configurations. This is an objective we will 
pursue in future work that we expect to report soon else- 
where. 
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